In [138]:
print(R.version)
file.path(R.home("bin"), "R")
In [139]:
library(IRdisplay) # display
library(fUnitRoots) # adfTest
library(TSA) # eacf
In [109]:
# detach("package:AFTSCode", unload = TRUE) # First detach the package
# unloadNamespace("AFTSCode")
In [140]:
library(AFTSCode)
In [ ]:
# plot_auto_arima_forecast_fig <- function(
# da_ts, eotr, h, npts, frequency,
# xreg=NULL, main=NULL, xlab=NULL, ylab=NULL, ylim = NULL, ...
# ) {
# par(bg = "white")
# # arima model
# if (is.null(seasonal)) {seasonal = list(order = c(0L, 0L, 0L), period = NA)} # default value
# tr_da_ts <- ts(da_ts[1:eotr], frequency = frequency, start = start(da_ts))
# if (is.null(xreg)) {
# tr_xreg <- NULL
# fc_xreg <- NULL
# } else {
# stopifnot("xreg should be of type matrix"=(is.matrix(xreg)))
# stopifnot("length(da_ts)!=dim(xreg)[1]"=(length(da_ts)==dim(xreg)[1]))
# tr_xreg <- xreg[1:eotr]
# fc_xreg <- xreg[(eotr+1):dim(xreg)[1]]
# }
# ts_fm <- auto.arima(tr_da_ts, xreg = tr_xreg, ...)
# print(ts_fm$nobs)
# # Forecast
# ts_fm$x <- tr_da_ts # https://stackoverflow.com/a/42464130/4307919
# ts_fc_res <- forecast(ts_fm, h = h, xreg = fc_xreg)
# # Plot forecast
# if (is.null(npts)) {
# npts <- eotr
# }
# xmin <- time(da_ts)[eotr] - npts / frequency
# xmax <- time(da_ts)[eotr] + (max(h, length(da_ts) - eotr) + 1) / frequency
# cat(xmin, ";", xmax)
# # # Label 1: Actual Observation Line
# plot(ts_fc_res, xlim = c(xmin, xmax), ylim = ylim, main = main, xlab = xlab, ylab = ylab)
# # Plot forecast mean (prepend the last observed data in the training dataset)
# dummy_1st_fmean_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$mean)), frequency = frequency, start = end(tr_da_ts))
# # # Label -: NULL
# lines(dummy_1st_fmean_ts)
# # # Label 2: Forecast Mean
# points(dummy_1st_fmean_ts, pch = 1)
# # Plot confidence interval (95%)
# dummy_1st_flower_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$lower[, 2])), frequency = frequency, start = end(tr_da_ts))
# dummy_1st_fupper_ts <- ts(c(c(da_ts[eotr]), as.numeric(ts_fc_res$upper[, 2])), frequency = frequency, start = end(tr_da_ts))
# # # Label 3: Forecast 95% Lower Bound
# lines(dummy_1st_flower_ts, lty = 2)
# # # Label 4: Forecast 95% Upper Bound
# lines(dummy_1st_fupper_ts, lty = 2)
# # Plot original data
# orig_plot_ts <- ts(da_ts[(eotr - npts + 1):length(da_ts)], frequency = frequency, start = time(da_ts)[eotr] - (npts - 1) / frequency)
# # # Label -: NULL
# lines(orig_plot_ts)
# # # Label 5: Actual Observation Points
# points(orig_plot_ts, pch = 19)
# legend(
# "topleft",
# legend = c(
# "Actual Obs Line", NULL, "Forecast Mean", "Forecast 95% Lower Bound",
# "Forecast 95% Upper Bound", NULL, "Actual Obs"
# ),
# # col = c("black", "red", "blue"),
# lty = c(1, NA, 2, 2, NA),
# pch = c(NA, 1, NA, NA, 19)
# )
# ts_fc_res
# }
2-1¶
- From Sec. 2.5.4
- $E[R_{101}|F_{100}] = c_0 - \theta_1 a_{100} = 0.2*0.01 = 0.002$
- $Var[R_{101}|F_{100}] = \sigma_{a}^2 = 0.025^2 = 0.000625$
- $E[R_{102}|F_{100}] = c_0 = 0$
- $Var[R_{102}|F_{100}] = (1+\theta_1^2)\sigma_{a}^2 = (1+0.2^2)*0.025^2 = 0.00065$
- From Sec. 2.5.1
- $\rho_1 = -(-0.2)/(1+0.2^2) = 0.192307692307692$
- $\rho_2 = 0$
2-2¶
- From Sec. 2.4.1
- $E[r_t]=\frac{\phi_0}{1-\phi_1}=0.01/(1-0.2)=0.0125$
- $Var[r_t]=\frac{\sigma_a^2}{1-\phi_1^2}=0.02/(1-0.2^2)=0.0208333333333333$
- From Sec. 2.4.2
- $\rho_0=1$
- $\rho_1=\phi_1\rho_0=0.2$
- $\rho_2=\phi_1\rho_1=0.04$
- From Sec. 2.4.4
- $E[r_{101}|F_{100}]=\phi_0+\phi_1*r_{100}=0.01+0.2*(-0.01)=0.008$
- $Var[r_{101}|F_{100}]=\sigma_a^2=0.02$
- $E[r_{102}|F_{100}]=\phi_0+\phi_1E[r_{101}|F_{100}]=0.01+0.2*0.008=0.0116$
- $Var[r_{102}|F_{100}]=Var[\phi_1 r_{101}+a_{102}|F_{100}]=Var[\phi_1 (\phi_0+\phi_1 r_{100}+a_{101})+a_{102}|F_{100}]=(1+\phi_1^2)\sigma_a^2=(1+0.2^2)*0.02=0.0208$
2-3¶
In [10]:
da = read.table("../AFTS_sol/data/m-unrate.txt", header = T)
da[1:5,]
In [15]:
unem_rate = da$Rate
lg_unem_rate = log(unem_rate)
unem_rate_ts = ts(unem_rate, frequency = 12, start = c(1948,1))
lg_unem_rate_ts = ts(lg_unem_rate, frequency = 12, start = c(1948,1))
In [143]:
plot_time_fig(unem_rate_ts, main = "Unemployment Rate", xlab = "Year")
pacf(unem_rate)
In [22]:
plot_time_fig(lg_unem_rate_ts, main = "log(Unemployment Rate)", xlab = "Year")
pacf(lg_unem_rate)
Unit-root test¶
We can reject the null-hypothesis from ADF test. In other words, the process doesn't have unit-root.- Seems like we cannot reject he null hypothesis when we use large enough lags in the ADF testing.
- Also from the ACF plots of the unemployment rate before and after taking difference, we can see that we need to take difference to cancel the strong serial correlation.
- From the ACF plot of the unemployment rate after taking difference, we see some damping sine and cosine wave, indicating that there could be some business cycle (pp42).
In [34]:
# "lags=12" is chosen from PACF plot
adf_test_res = adfTest(unem_rate_ts, lags = 12, type = c("ct"))
adf_test_res
In [31]:
lag.max = 60
plot_acf(da = unem_rate, lag.max = lag.max, main = "Unemployment Rate")
pacf(unem_rate, lag.max = lag.max)
In [35]:
diff_unem_rate = diff(unem_rate)
diff_unem_rate_ts = ts(diff_unem_rate, frequency = 12, start = c(1948,2))
In [36]:
plot_acf(da = diff_unem_rate, lag.max = lag.max, main = "diff(Unemployment Rate)")
pacf(diff_unem_rate, lag.max = lag.max)
In [144]:
plot_acf(da = diff(unem_rate, 12), lag.max = lag.max, main = "diff(Unemployment Rate, 12)")
pacf(diff(unem_rate, 12), lag.max = lag.max)
Order determination¶
- Try to use ARIMA(4, 1, 5) model.
In [172]:
perform_and_print_eacf <- function(da, ar.max, ma.max) {
eacf_obj <- eacf(da, ar.max = ar.max, ma.max = ma.max)
eacf_stats_tb = format(as.data.frame(eacf_obj$eacf), digits = 3)
names(eacf_stats_tb) <- seq(from = 0, to = ma.max)
display(eacf_stats_tb)
display(eacf_obj$symbol)
# pp67, asymptotic standard error of EACF
display(2/sqrt(length(da)))
c(eacf_obj, eacf_stats_tb)
}
In [173]:
unem_rate_eacf_res <- perform_and_print_eacf(diff_unem_rate, ar.max = 25, ma.max = 12)
Fit an ARIMA(3,1,5) model¶
In [108]:
unem_rate_mod = stats::arima(unem_rate_ts, order = c(3,1,5), include.mean = F)
unem_rate_mod
In [110]:
unem_rate_mod_1 = arima(unem_rate_ts, order = c(3,1,5))
unem_rate_mod_1
In [111]:
diff_unem_rate_mod = arima(diff_unem_rate_ts, order = c(3,0,5), include.mean = F)
diff_unem_rate_mod
In [112]:
plot_time_fig(unem_rate_mod$residuals, main = "Residuals")
plot_acf(unem_rate_mod$residuals, main = "ACF of Residuals")
Box-Ljung test¶
- pp33 very briefly talks about how to choose lag in the
Box.test - From the following testing result, we cannot reject the null-hypothesis. So the model specified is adequate.
In [113]:
length(unem_rate); log(length(unem_rate))
In [114]:
Box.test(unem_rate_mod$residuals, lag = 7, type = 'Ljung')
In [115]:
Box.test(unem_rate_mod$residuals, lag = 12, type = 'Ljung')
- Diagnose time-series ARIMA model
In [161]:
help(tsdiag)
In [162]:
par(bg = 'white')
tsdiag(unem_rate_mod, gof.lag=36)
Check business cycle¶
- pp42
- Compute the average length of business cycles
- Doesn't matter if using the characteristic roots or the inverse of them.
- Unit is "month". Average business cycle is about 14.2 months, which seems to be shoter than expected from the time plot above.
- The business cycle is very sensitive to the order chosen in fitting ARIMA model. When I use ARIMA(2,0,5) or ARIMA(4,0,5), the business cycle doens't exist or too short.
In [69]:
find('arima')
In [79]:
help("arima", package = "TSA")
In [80]:
help("arima", package = "stats")
In [116]:
sqrt(unem_rate_mod$sigma2)
In [117]:
unem_rate_mod$coef
In [118]:
ar_poly = c(1, -unem_rate_mod$coef[1:3]) # Characteristic equation for AR process
roots = polyroot(ar_poly)
roots
In [119]:
for (i in 1:3) {
print(paste("====", i))
print(as.numeric(Mod(1-unem_rate_mod$coef[1]*roots[i]-unem_rate_mod$coef[2]*roots[i]^2-unem_rate_mod$coef[3]*roots[i]^3-unem_rate_mod$coef[4]*roots[i]^4)))
print(as.numeric(Mod(roots[i])))
}
- Modulus of roots
In [120]:
roots; Mod(roots)
In [121]:
1/roots; Mod(1/roots)
- Compute the average length of business cycles
- Doesn't matter if using the characteristic roots or the inverse of them.
In [122]:
roots[1]; Re(roots[1]); Im(roots[1])
In [125]:
k = 2*pi/acos(Re(roots[1])/Mod(roots[1]))
k
In [126]:
k = 2*pi/acos(Re(1/roots[1])/Mod(1/roots[1]))
k
Forecast¶
In [128]:
length(unem_rate_ts); unem_rate_ts
In [140]:
npts = 8
eotr = length(unem_rate_ts)
h = 4
freq = 12
order = c(3,1,5)
fixed = NULL
seasonal = list(order = c(0L, 0L, 0L), period = NA)
unem_rate_fc_res = plot_arima_forecast_fig(
da_ts=unem_rate_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=F, transform.pars=NULL,
main="Forecasts from ARIMA(3,1,5) for unemp_rate",
xlab="Year", ylab="Unemp-Rate", ylim=c(5, 11)
)
Alternative solutions (only use AR model)¶
In [147]:
par(bg = 'white')
# decay slowly
acf(unem_rate)
par(mfrow = c(2, 1))
acf(diff(unem_rate))
pacf(diff(unem_rate))
acf(diff(unem_rate, 12))
pacf(diff(unem_rate, 12))
cat("order =", ar(unem_rate, method = 'mle')$order, "\n")
m1 = arima(unem_rate_ts, order = c(11, 0, 0))
m1
# Follow the pp51 to modify and refine the model based on previous fit
m1 = arima(unem_rate_ts, order = c(11, 0, 0), fixed = c(NA, NA, 0, 0, 0, NA, 0, 0, 0, NA, NA, NA))
m1
# FIXME
m2 = arima(unem_rate_ts, order = c(2, 1, 1), seasonal = list(order = c(1, 0, 1), period = 12))
m2
tsdiag(m1, gof = 36)
tsdiag(m2, gof = 36)
predict(m1, 4)
predict(m2, 4)
In [149]:
ar_poly_1 = c(1, -m1$coef[1:11]) # Characteristic equation for AR process
roots_1 = polyroot(ar_poly_1)
roots_1
In [155]:
Im(roots_1[8])
In [156]:
for (i in 1:length(roots_1)) {
if (abs(Im(roots_1[i]))>1e-8) {
print(i)
print(roots_1[i])
print(2*pi/acos(Re(roots_1[i])/Mod(roots_1[i])))
}
}
In [151]:
ar_poly_2 = c(1, -m2$coef[1:2]) # Characteristic equation for AR process
roots_2 = polyroot(ar_poly_2)
roots_2
2-4¶
In [3]:
da = read.table("../AFTS_data/Ch02/m-deciles08.txt", header = T)
dim(da); da[1:5,]
In [4]:
dec2_ts = ts(da$CAP2RET, frequency = 12, start = c(1970, 1))
dec10_ts = ts(da$CAP10RET, frequency = 12, start = c(1970, 1))
In [5]:
plot_time_fig(dec2_ts, main = "CRSP Decile Index", xlab = "Year")
lines(dec10_ts, type = , col = 'red', lty = 2)
legend(
"topright",
legend = c("Decile 2", "Decile 10"),
col = c("black", "red"),
lty = c(1, 2),
pch = c(8, NA)
)
(a) Test autocorrelations¶
- Reject the null-hypothesis (i.e., the first 12 lags have no serial correlation) for Decile 2. So for Decile 2, the first 12 lags have serial correlation.
- Cannot reject the null-hypothesis for Decile 10. So for Decile 10, the first 12 lags have no significant serial correlation.
In [6]:
par(mfrow = c(2, 1), bg = 'white')
acf(da$CAP2RET)
pacf(da$CAP2RET)
Box.test(da$CAP2RET, lag = 12, type = 'Ljung')
In [7]:
par(mfrow = c(2, 1), bg = 'white')
acf(da$CAP10RET)
pacf(da$CAP10RET)
Box.test(da$CAP10RET, lag = 12, type = 'Ljung')
(b) Fit a model and perform model checking¶
In [8]:
dec2_eacf_res <- perform_and_print_eacf(da$CAP2RET, ar.max = 25, ma.max = 12)
- Fit
ARIMA(1,0,1)model, b/c only this model has all coefs significant.
In [10]:
dec2_mod_1 <- arima(da$CAP2RET, order = c(1, 0, 1), include.mean = T)
dec2_mod_1; Box.test(dec2_mod_1$residuals, lag = 12, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_1, gof.lag = 36)
In [11]:
# Follow the pp51 to modify and refine the model based on previous fit
dec2_mod_2 <- arima(da$CAP2RET, order = c(0, 0, 1), include.mean = T)
dec2_mod_2; Box.test(dec2_mod_2$residuals, lag = 10, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_2, gof.lag = 36)
In [12]:
dec2_mod_3 <- arima(da$CAP2RET, order = c(4, 0, 3), include.mean = T)
dec2_mod_3; Box.test(dec2_mod_3$residuals, lag = 12, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_3, gof.lag = 36)
In [13]:
dec2_mod_3 <- arima(
da$CAP2RET,
order = c(4, 0, 3),
# Follow the pp51 to modify and refine the model based on previous fit
fixed = c(0, NA, 0, NA, 0, NA, 0, NA),
include.mean = T
)
dec2_mod_3; Box.test(dec2_mod_3$residuals, lag = 12, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_3, gof.lag = 36)
(c) Forecast using ARMA model¶
In [22]:
npts = 12
eotr = length(dec2_ts)-npts
h = npts
freq = 12
order = c(0,0,1)
fixed = NULL
seasonal = NULL
dec2_fc_res = plot_arima_forecast_fig(
da_ts=dec2_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=T, transform.pars=NULL,
main="Forecasts from ARIMA(0,0,1) for Decile 2 index",
xlab="Year", ylab="Unemp-Rate"# , ylim=c(5, 11)
)
dec2_fc_tb = comb_forecast_res(dec2_fc_res, dec2_ts, eotr, freq)
dec2_fc_tb
In [17]:
class(dec2_fc_res); attributes(dec2_fc_res); dec2_fc_res$model
In [20]:
help("predict")
In [21]:
predict(dec2_fc_res$model, 12)
In [23]:
dec2_fc_tb
2-5¶
- There is some long-range dependency.
In [25]:
da = read.table("../AFTS_sol/data/d-ibm3dx7008.txt", header = T)
ibm = da$rtn
plot_acf(abs(ibm), lag.max = 100)
2-6¶
- This is a multiplicative seasonal model.
- The seasonal prediction is very good.
- The data set in this problem bears a lot similarity with the JnJ quarterly earning data set in Sec. 2.8. After taking one regular difference, the seasonality becomes very obvious. This means that after taking one regular difference, we need to model the AR or MA pattern of the resulting time series.
- The logic of determining the orders in arima model from the various
pacfandacfplots can be as following.- The plots for
powshow strong serial correlation. So we need to take diff. - The plots for
diff(pow)show strong seasonality in theacf. So we need to take the seasonal diff. - The plots for
diff(pow, 12)show some autocorrelation in theacf. So we need to take both diff and seasonal diff. - The plots for
diff(diff(pow), 12)show some autocorrelation in theacf. So we need to apply MA model ondiff(diff(pow), 12). Thus this is a multiplicative seasonal model.
- The plots for
In [100]:
da = read.table("../AFTS_sol/data/power6.txt", header = F)
dim(da); da[1:5,]
In [101]:
par(bg = 'white')
pow = da[, 1]
pow_ts = ts(pow, frequency = 12, start = c(2000, 1))
# From the time plot, should immediately think that we can apply multiplicative seasonal model
plot(pow_ts, type = 'o', xlab = 'Year', ylab = 'Power', main = "Power Consumption")
plot(diff(pow_ts), type = 'o', xlab = 'Year', ylab = 'diff(Power)', main = "Power Consumption")
plot(diff(pow_ts, 12), type = 'o', xlab = 'Year', ylab = 'diff(Power, 12)', main = "Power Consumption")
plot(diff(diff(pow_ts), 12), type = 'o', xlab = 'Year', ylab = 'diff(diff(Power), 12)', main = "Power Consumption")
In [102]:
par(mfrow = c(2, 1), bg = 'white')
pacf(pow)
acf(pow)
pacf(diff(pow))
acf(diff(pow))
pacf(diff(pow, 12))
acf(diff(pow, 12))
pacf(diff(diff(pow), 12))
acf(diff(diff(pow), 12))
In [103]:
# The eacf is not helpful in determining the order of ARIMA model
pow_eacf_res <- perform_and_print_eacf(pow, ar.max = 25, ma.max = 12)
In [104]:
par(bg = 'white')
# FIXME
m1 = arima(pow, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
m1
tsdiag(m1, gof = 36)
# predict(m1, 24)
In [111]:
pow_ts = ts(pow, frequency = 12, start = c(2010, 1))
npts = 24
eotr = length(pow_ts)-npts
h = npts
freq = 12
# The order and seasonal following in combination gives a multiplicative seasonal model.
order = c(0, 1, 1)
seasonal = list(order = c(0, 1, 1), period = 12)
fixed = NULL
pow_fc_res = plot_arima_forecast_fig(
da_ts=pow_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=T, transform.pars=NULL,
main="Forecasts from ARIMA(0,1,1)(0,1,1)[12] for Power consumption",
xlab="Year", ylab="Power", ylim=c(9.6, 10.1)
)
In [44]:
pow_fc_tb = comb_forecast_res(pow_fc_res, pow_ts, eotr, freq)
pow_fc_tb
2-7 (Equal-Weighted Index)¶
In [114]:
require(xts)
require(sandwich)
require(lmtest)
require(forecast)
In [13]:
find("read.table")
Fit a time series model for CRSP equal-weighted index¶
- Test the effect of trading days using regression model without time series error.
- Add time series error model by using ARIMA model. The order of ARIMA can be seen from some acf and pacf plots; Also can look at eacf plots and tables. After getting a starting point for the
orderparameter ofarimafunction, we can trim the number of coefs by reducing element numbers inorder, until we have all coefs close to be statistically significant. - The logic of determining the order in
arimamodel is as following.- We see that the
m1$residualshave some autocorrelation inacf. We need to take diff. - After taking one diff, we see that the
diff(m1$residuals)have some autocorrelation in bothpacfandacf. We take seasonal diff. - After taking one seasonal diff, we see that
pacfhas some strong seasonality. Andacfhas some autocorrelation. Let's look at effect of both regular diff and seasonal diff. - After taking both regular and seasonal diff, we see that the autocorrelation in both
pacfandacfgetting stronger. We probably don't want to do this. - So, we need to apply seasonal diff and model the decaying autocorrelation in the
pacf. That means modeling the AR process.
- We see that the
In [113]:
da = read.table("../AFTS_sol/data/d-ibm3dxwkdays8008.txt", header = T)
da[1:5,]
In [115]:
ew = da$ew * 100
ew_ts = ts(ew, frequency = 252, start = c(1980, 1, 2))
plot(
xts(ew, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'ew'
)
- Trading day effects are all significant.
In [116]:
M = da$M
T = da$T
W = da$W
R = da$R
m1 = lm(ew ~ M + T + W + R)
summary(m1)
- Trading day effects are still significant. So HAC estimator doesn't change the conclusion of weekday effects.
In [17]:
find("NeweyWest")
In [117]:
coeftest(m1, NeweyWest(m1, lag = 12, prewhite = F))
- Residuals have serial correlation.
- Box-Ljung test tells us the null-hypothesis can be rejected.
In [118]:
Box.test(m1$residuals, lag = 12, type = 'Ljung')
- Residuals have some autocorrelations.
In [119]:
par(bg = 'white')
plot_time_fig(m1$residuals)
plot_time_fig(m1$residuals[1:200])
plot_time_fig(m1$residuals[(length(m1$residuals)-200):length(m1$residuals)])
In [122]:
par(mfrow = c(2, 1), bg = 'white')
pacf(m1$residuals)
acf(m1$residuals)
pacf(diff(m1$residuals))
acf(diff(m1$residuals))
pacf(diff(m1$residuals, 5))
acf(diff(m1$residuals, 5))
pacf(diff(diff(m1$residuals), 5))
acf(diff(diff(m1$residuals), 5))
- Fit an ARIMA model with time series error.
In [59]:
ew_eacf_res <- perform_and_print_eacf(ew, ar.max = 25, ma.max = 12)
In [56]:
par(bg = 'white')
m2 = arima(
ew, order = c(2, 0, 2),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
m2
tsdiag(m2, gof = 20)
In [58]:
par(bg = 'white')
m3 = arima(
ew, order = c(2, 0, 3),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
m3
tsdiag(m3, gof = 20)
In [24]:
par(bg = 'white')
m4 = arima(
ew, order = c(4, 0, 4),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
m4
tsdiag(m4, gof = 20)
In [25]:
par(bg = 'white')
m4 = arima(
ew, order = c(3, 0, 3),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
m4
tsdiag(m4, gof = 20)
In [126]:
par(bg = 'white')
m5 = arima(
ew, order = c(2, 0, 2),
# seasonal = list(order = c(0, 0, 1), period = 5),
xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
In [125]:
par(bg = 'white')
m5 = arima(
ew, order = c(2, 0, 2),
seasonal = list(order = c(0, 0, 1), period = 5),
xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
In [128]:
par(bg = 'white')
m5 = arima(
ew, order = c(2, 0, 2),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
In [127]:
par(bg = 'white')
m5 = arima(
ew, order = c(2, 0, 2),
seasonal = list(order = c(1, 0, 1), period = 5),
xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
Forecast using the model¶
In [7]:
help(auto.arima)
Test auto.arima with a subset of time series data¶
In [17]:
da[6560:6570,]
In [19]:
sub_ew_ts = ts(ew[6565:length(ew)], frequency = 252, start = c(2006, 1, 3))
sub_da = da[6565:length(ew),]
length(sub_ew_ts)
In [23]:
sub_ts_fm <- auto.arima(
sub_ew_ts,
d = 0,
D = 0,
max.p = 2,
max.q = 2,
max.P = 1,
max.Q = 0,
max.order = 5,
# nmodels = 10,
seasonal = TRUE,
method = "ML",
allowmean = TRUE,
xreg = as.matrix(sub_da[, 8:11]),
# stepwise = FALSE,
# parallel = TRUE,
# num.cores = 6
)
sub_ts_fm
In [25]:
class(sub_ts_fm); attributes(sub_ts_fm)
In [29]:
h = 20
sub_ts_fc_res <- forecast(sub_ts_fm, h = h, xreg = as.matrix(sub_da[(dim(sub_da)[1]-h):(dim(sub_da)[1]), 8:11]))
summary(sub_ts_fc_res)
In [31]:
par(bg = 'white')
plot(sub_ts_fc_res)
In [32]:
sub_ts_fc_res$mean
In [33]:
sub_ts_fc_res$lower; sub_ts_fc_res$upper
Plot forecast results from auto.arima model with a subset of time series¶
In [41]:
help(stopifnot)
In [112]:
npts = 20
eotr = length(sub_ew_ts)-npts
h = npts
freq = 252
xreg = as.matrix(sub_da[, 8:11])
sub_ew_fc_res = plot_auto_arima_forecast_fig(
da_ts=sub_ew_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
xreg=xreg,
main="Forecasts from ARIMA(1,1,0)\nfor CRSP Equal-Weighted Index",
xlab="Year", ylab="CRSP Equal-Weighted Index"# , ylim=c(5, 11)
)
In [49]:
sub_ew_fc_tb = comb_forecast_res(sub_ew_fc_res, sub_ew_ts, eotr, freq)
sub_ew_fc_tb
In [62]:
class(sub_ew_fc_res); attributes(sub_ew_fc_res);
class(sub_ew_fc_res$model); attributes(sub_ew_fc_res$model);
In [59]:
sub_ew_fc_res$model
All data takes too long time to run auto.arima¶
Dummy example for auto.arima¶
In [50]:
# Assuming 'forecast' package is installed and loaded
library(forecast)
# Simulate some time series data
set.seed(123)
ts_data <- ts(rnorm(100), frequency = 12, start = c(2000, 1))
# Simulate exogenous regressor data for model fitting
xreg_train <- matrix(rnorm(200), ncol = 2) # Two exogenous regressors
# Fit an ARIMA model with exogenous regressors
# fit <- auto.arima(ts_data, xreg = xreg_train)
fit <- auto.arima(ts_data, xreg = xreg_train)
# Simulate future values of the exogenous regressors for forecasting
xreg_future <- matrix(rnorm(24), ncol = 2) # Two exogenous regressors for 12 future periods
# Forecast using the model and future values of exogenous regressors
forecast_result <- forecast(fit, xreg = xreg_future)
# Plot the forecast
par(bg = 'white')
plot(forecast_result)
In [57]:
npts = 24
eotr = length(ts_data)-npts
h = npts
freq = 12
xreg = xreg_train
sim_data_fc_res = plot_auto_arima_forecast_fig(
da_ts=ts_data, eotr=eotr, h=h, npts=npts, frequency=freq,
xreg=xreg,
main="Forecasts from ARIMA(0,0,0) for Simulated Data",
xlab="Year", ylab="Simulated Data"# , ylim=c(5, 11)
)
2-8 (S&P)¶
In [141]:
require(xts)
In [145]:
da = read.table("../AFTS_sol/data/d-ibm3dxwkdays8008.txt", header = TRUE)
dim(da); da[1:5,]
In [155]:
plot_pacf_acf(sp, lag.max = 20, main = "SP")
plot_pacf_acf(diff(sp), lag.max = 20, main = "diff(SP)")
plot_pacf_acf(diff(sp, 5), lag.max = 20, main = "diff(SP, 5)")
plot_pacf_acf(diff(diff(sp), 5), lag.max = 20, main = "diff(diff(SP), 5)")
(a)¶
- No Friday effect.
In [159]:
sp = da$sp * 100
sp_ts = ts(sp, frequency = 252, start = c(1980, 1, 2))
plot(
xts(sp, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'sp'
)
M = da$M
T = da$T
W = da$W
R = da$R
F = da$F
In [162]:
sp_m1 = lm(sp ~ M + T + W + R + F + 0)
summary(sp_m1);
In [163]:
Box.test(sp_m1$residuals, lag = 12, type = 'Ljung');
In [220]:
sp_m1 = lm(sp ~ 1 + T + W)
summary(sp_m1);
(b)¶
- There are some serial correlation in the residuals.
In [164]:
plot_pacf_acf(sp_m1$residuals, lag.max = 20, main = 'obs')
plot_pacf_acf(diff(sp_m1$residuals), lag.max = 20, main = 'diff(obs)')
plot_pacf_acf(diff(sp_m1$residuals, 5), lag.max = 20, main = "diff(obs, 5)")
plot_pacf_acf(diff(diff(sp_m1$residuals), 5), lag.max = 20, main = "diff(diff(obs), 5)")
In [169]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(0, 0, 1),
# seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [170]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(1, 0, 0),
# seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [171]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(2, 0, 0),
# seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [172]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(2, 0, 0),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [174]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(3, 0, 1),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [195]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(4, 0, 0),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [196]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(4, 0, 0),
fixed = c(NA, NA, 0, NA, 0, 0, 0, 0, 0),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [198]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(5, 0, 0),
fixed = c(NA, NA, 0, NA, 0, 0, 0, 0, 0, 0),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [ ]:
In [200]:
par(bg = 'white')
sp_m3 = arima(
sp, order = c(2, 0, 0),
xreg = da[, 8:11]
)
sp_m3
tsdiag(sp_m3, gof = 20)
In [201]:
par(bg = 'white')
sp_m3 = arima(
sp, order = c(2, 0, 2),
xreg = da[, 8:11]
)
sp_m3
tsdiag(sp_m3, gof = 20)
In [203]:
par(bg = 'white')
sp_m3 = arima(
sp, order = c(2, 0, 2),
fixed = c(0, NA, 0, NA, 0, 0, 0, 0, 0),
xreg = da[, 8:11]
)
sp_m3
tsdiag(sp_m3, gof = 20)
In [ ]:
In [ ]:
In [213]:
sp_m4 = arima(sp, order = c(2, 0, 2), xreg = da[, 8:12], include.mean = FALSE)
sp_m4;
par(bg = 'white')
tsdiag(sp_m4, gof = 20)
- Adding T and W, b/c they show significance in linear model.
In [215]:
sp_m4 = arima(sp, order = c(2, 0, 2), xreg = da[, c(9, 10, 12)], include.mean = FALSE,
fixed = c(0, NA, 0, NA, NA, NA, NA))
sp_m4;
par(bg = 'white')
tsdiag(sp_m4, gof = 20)
Box.test(sp_m4$residuals, lag = 10, type = 'Ljung')
In [ ]:
2-9 (IBM)¶
In [217]:
require(xts)
In [219]:
da = read.table("../AFTS_sol/data/d-ibm3dxwkdays8008.txt", header = TRUE)
da[1:5,]
(a)¶
- No Friday effect.
In [221]:
ibm = da$ibm * 100
plot(xts(ibm, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'ibm')
M = da$M
T = da$T
W = da$W
R = da$R
F = da$F
In [223]:
ibm_m1 = lm(ibm ~ M + T + W + R + F + 0)
summary(ibm_m1)
In [225]:
Box.test(ibm_m1$residuals, lag = 12, type = 'Ljung')
(b)¶
- There is no much serial correlation in the residuals.
In [231]:
lag.max = 20
par(mfrow = c(2, 1), bg = 'white')
pacf(ibm_m1$residuals, lag.max =lag.max)
acf(ibm_m1$residuals, lag.max =lag.max)
pacf(diff(ibm_m1$residuals), lag.max =lag.max)
acf(diff(ibm_m1$residuals), lag.max =lag.max)
pacf(diff(ibm_m1$residuals, 5), lag.max =lag.max)
acf(diff(ibm_m1$residuals, 5), lag.max =lag.max)
pacf(diff(diff(ibm_m1$residuals), 5), lag.max =lag.max)
acf(diff(diff(ibm_m1$residuals), 5), lag.max =lag.max)
(c)¶
In [228]:
ibm_m2 = arima(ibm, order = c(0, 0, 1), xreg = da[, 8:12], include.mean = FALSE)
ibm_m2;
par(bg = 'white')
tsdiag(ibm_m2, gof = 35)
In [229]:
Box.test(ibm_m2$residuals, lag = 12, type = 'Ljung')
In [234]:
ibm_m2 = arima(ibm, order = c(0, 0, 1), xreg = da[, c(8, 9)], include.mean = TRUE)
ibm_m2;
par(bg = 'white')
tsdiag(ibm_m2, gof = 35)
In [233]:
ibm_m2 = arima(ibm, order = c(0, 0, 1), xreg = da[, c(8, 9)], include.mean = FALSE)
ibm_m2;
par(bg = 'white')
tsdiag(ibm_m2, gof = 35)
2-10¶
- pp10: Test statistics for skewness and kurtosis.
- They are all significant.
In [237]:
require(fBasics)
In [235]:
da1 = read.table("../AFTS_sol/data/w-Aaa.txt", header = FALSE)
da2 = read.table("../AFTS_sol/data/w-Baa.txt", header = FALSE)
da1[1:5,]; da2[1:5,];
In [238]:
Aaa = da1[, 4]
Baa = da2[, 4]
apply(cbind(Aaa, Baa), 2, basicStats)
ts = skewness(Aaa) / sqrt(6 / length(Aaa))
cat("ts =", ts, "\n")
ps = (1 - pnorm(ts)) * 2
cat("ps =", ps, "\n")
tk = kurtosis(Aaa) / sqrt(24 / length(Aaa))
cat("tk =", tk, "\n")
pk = (1 - pnorm(tk)) * 2
cat("pk =", pk, "\n")
ts = skewness(Baa) / sqrt(6 / length(Baa))
cat("ts =", ts, "\n")
ps = (1 - pnorm(ts)) * 2
cat("ps =", ps, "\n")
tk = kurtosis(Baa) / sqrt(24 / length(Baa))
cat("tk =", tk, "\n")
pk = (1 - pnorm(tk)) * 2
cat("pk =", pk, "\n")
2-11¶
In [239]:
da1 = read.table("../AFTS_sol/data/w-Aaa.txt", header = FALSE)
da1[1:5,];
In [252]:
aaa = da1[,4]
aaa_ts = ts(aaa, frequency = 52, start = c(1962, 1))
lg_aaa = log(aaa)
lg_aaa_ts = ts(lg_aaa, frequency = 52, start = c(1962, 1))
Fitting solution¶
In [253]:
sub_sz = 100
plot_time_fig(lg_aaa_ts, main = "AAA Bond", xlab = "Year", ylab = "lg(Bond Weekly Yield)")
plot_time_fig(lg_aaa_ts[1:sub_sz], main = "AAA Bond", xlab = "Year", ylab = "lg(Bond Weekly Yield)")
plot_time_fig(lg_aaa_ts[(length(lg_aaa_ts)-sub_sz):length(lg_aaa_ts)], main = "AAA Bond", xlab = "Year", ylab = "lg(Bond Weekly Yield)")
In [243]:
sub_sz = 100
plot_time_fig(aaa_ts, main = "AAA Bond", xlab = "Year", ylab = "Bond Weekly Yield")
plot_time_fig(aaa_ts[1:sub_sz], main = "AAA Bond", xlab = "Year", ylab = "Bond Weekly Yield")
plot_time_fig(aaa_ts[(length(aaa_ts)-sub_sz):length(aaa_ts)], main = "AAA Bond", xlab = "Year", ylab = "Bond Weekly Yield")
In [250]:
lag.max = 52
par(mfrow = c(2, 1), bg = 'white')
pacf(aaa, lag.max = lag.max)
acf(aaa, lag.max = lag.max)
pacf(diff(aaa), lag.max = lag.max)
acf(diff(aaa), lag.max = lag.max)
pacf(diff(aaa, 4), lag.max = lag.max)
acf(diff(aaa, 4), lag.max = lag.max)
pacf(diff(aaa, 13), lag.max = lag.max)
acf(diff(aaa, 13), lag.max = lag.max)
pacf(diff(aaa, 26), lag.max = lag.max)
acf(diff(aaa, 26), lag.max = lag.max)
pacf(diff(diff(aaa), 26), lag.max = lag.max)
acf(diff(diff(aaa), 26), lag.max = lag.max)
In [ ]:
In [256]:
par(bg = 'white')
aaa_m1 = arima(
aaa, order = c(0, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m1
tsdiag(aaa_m1, gof = 52)
In [259]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(1, 1, 0),
seasonal = list(order = c(1, 1, 0), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [258]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(1, 1, 1),
seasonal = list(order = c(1, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [260]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(1, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [261]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(2, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [262]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(3, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [263]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(4, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [265]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(9, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [266]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(9, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26),
fixed = c(NA, 0, NA, 0, NA, NA, NA, NA, NA, 0, NA)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [267]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(9, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26),
fixed = c(NA, 0, NA, 0, NA, 0, NA, 0, NA, 0, NA)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [273]:
par(bg = 'white')
aaa_m2 = arima(aaa, order = c(9, 1, 1))
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [276]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(9, 1, 0),
fixed = c(NA, 0, NA, 0, 0, 0, NA, 0, NA)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [ ]:
Alternative solution¶
In [268]:
plot(
xts(aaa, order.by = as.Date(paste(da1[, 1], da1[, 2], da1[, 3], sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'Aaa'
)
In [270]:
cat("order =", ar(diff(aaa), method = 'mle')$order, "\n")
In [271]:
aaa_al_m1 = arima(aaa, order = c(9, 1, 0))
aaa_al_m1;
par(bg = 'white')
tsdiag(aaa_al_m1, gof = 52)
In [291]:
aaa_al_m1 = arima(
aaa, order = c(9, 1, 0),
seasonal = list(order = c(1, 0, 0), period = 13)
)
aaa_al_m1;
par(bg = 'white')
tsdiag(aaa_al_m1, gof = 52)
In [277]:
aaa_al_m2 = arima(aaa, order = c(2, 1, 3))
aaa_al_m2;
par(bg = 'white')
tsdiag(aaa_al_m2, gof = 52)
In [288]:
aaa_al_m2 = arima(
aaa, order = c(2, 1, 3),
seasonal = list(order = c(1, 0, 0), period = 26)
)
aaa_al_m2;
par(bg = 'white')
tsdiag(aaa_al_m2, gof = 52)
Forecast¶
In [282]:
npts = 26
eotr = length(aaa_ts)-npts
h = npts
freq = 52
order = c(9,1,0)
fixed = NULL
seasonal = NULL
aaa_fc_res = plot_arima_forecast_fig(
da_ts=aaa_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=TRUE, transform.pars=NULL,
main="Forecasts from ARIMA(9,1,0) for Aaa weekly bond yields",
xlab="Year", ylab="Aaa weekly bond yields", ylim=c(4, 8)
)
In [283]:
aaa_fc_tb = comb_forecast_res(aaa_fc_res, aaa_ts, eotr, freq)
aaa_fc_tb
In [293]:
npts = 26
eotr = length(aaa_ts)-npts
h = npts
freq = 52
order = c(2,1,3)
fixed = NULL
seasonal = NULL
aaa_fc_res_1 = plot_arima_forecast_fig(
da_ts=aaa_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=TRUE, transform.pars=NULL,
main="Forecasts from ARIMA(2,1,3) for Aaa weekly bond yields",
xlab="Year", ylab="Aaa weekly bond yields", ylim=c(4, 8)
)
In [285]:
aaa_fc_tb_1 = comb_forecast_res(aaa_fc_res_1, aaa_ts, eotr, freq)
aaa_fc_tb_1
In [298]:
npts = 26
eotr = length(aaa_ts)-npts
h = npts
freq = 52
order = c(2,1,3)
fixed = NULL
seasonal = list(order = c(0, 0, 1), period = 13)
aaa_fc_res_2 = plot_arima_forecast_fig(
da_ts=aaa_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=TRUE, transform.pars=NULL,
main="Forecasts from ARIMA(2,1,3)(0,0,1)[13] for Aaa weekly bond yields",
xlab="Year", ylab="Aaa weekly bond yields", ylim=c(4, 8)
)
In [ ]:
2-12¶
In [299]:
da1 = read.table("../AFTS_sol/data/w-Aaa.txt", header = FALSE)
da2 = read.table("../AFTS_sol/data/w-Baa.txt", header = FALSE)
da1[1:5,]; da2[1:5,];
In [302]:
freq = 52
start = c(1962, 1)
aaa = da1[,4]
aaa_ts = ts(aaa, frequency = freq, start = start)
bbb = da2[,4]
bbb_ts = ts(bbb, frequency = freq, start = start)
aaac = diff(aaa)
bbbc = diff(bbb)
aaac_ts = ts(aaac, frequency = freq, start = start)
bbbc_ts = ts(bbbc, frequency = freq, start = start)
In [316]:
apply(cbind(aaac, bbbc), 2, basicStats)
Exploration¶
In [301]:
par(bg = 'white')
plot(aaa_ts, main = "Weekly Aaa or Bbb Bond Yield", xlab = "Year", ylab = "Yield")
lines(bbb_ts, lty = 2)
legend(
"topleft",
legend = c("Aaa", "Bbb"),
lty = c(1, 2)
)
In [304]:
par(bg = 'white')
plot(aaac_ts, main = "Weekly Aaa Bond Yield", xlab = "Year", ylab = "Yield")
plot(bbbc_ts, main = "Weekly Bbb Bond Yield", xlab = "Year", ylab = "Yield")
In [305]:
par(bg = "white")
plot(x = aaa, y = bbb, main = "Aaa vs Bbb", xlab = "Aaa Yield", ylab = "Bbb Yield")
plot(x = aaac, y = bbbc, main = "Aaa Change vs Bbb Change", xlab = "diff(Aaa Yield)", ylab = "diff(Bbb Yield)")
Fitting¶
- Strong serial correlation in ACF of residuals, indicating that there is a unit-root process.
- Also conducted ADF unit-root tests. The results show that null-hypothesis cannot be rejected.
- Null-hypothesis is that the time series has unit-root.
- Box-Ljung test shows that numm-hypothesis can be rejected.
- Null-hypothesis is that the time series at certain lag has no autocorrelation.
- After taking one diff, the autocorrelation is largely gone. But Box-Ljung test is still significant. We can model the residuals as a time series.
Fitting the bond yields¶
In [306]:
bond_m1 = lm(aaa~bbb)
summary(bond_m1)
In [308]:
par(bg = 'white')
resi_ts = ts(bond_m1$residuals, frequency = freq, start = start)
plot(resi_ts, type = 'l', main = "Simple Linear Model", xlab = "Year")
lag.max = 36
par(mfrow = c(2,1))
pacf(bond_m1$residuals, lag.max = lag.max, main = "Residuals PACF Plot")
acf(bond_m1$residuals, lag.max = lag.max, main = "Residuals ACF Plot")
In [309]:
adf_test_res_a = adfTest(aaa, lags = 12, type = c("ct"))
adf_test_res_b = adfTest(bbb, lags = 12, type = c("ct")) # 2.7.3. Trend-Stationary Time-Series
adf_test_res_a; adf_test_res_b
In [310]:
Box.test(bond_m1$residuals, lag = 12, type = 'Ljung')
Fitting the bond yield changes¶
In [311]:
bond_m2 = lm(aaac~-1+bbbc) # No intercept
summary(bond_m2)
In [312]:
par(bg = 'white')
resi_ts = ts(bond_m2$residuals, frequency = freq, start = start)
plot(resi_ts, type = 'l', main = "Simple Linear Model", xlab = "Year")
lag.max = 36
par(mfrow = c(2,1))
pacf(bond_m2$residuals, lag.max = lag.max, main = "Residuals PACF Plot")
acf(bond_m2$residuals, lag.max = lag.max, main = "Residuals ACF Plot")
In [314]:
adf_test_res_ac = adfTest(aaac, lags = 12, type = c("ct"))
adf_test_res_bc = adfTest(bbbc, lags = 12, type = c("ct")) # 2.7.3. Trend-Stationary Time-Series
adf_test_res_ac; adf_test_res_bc
In [313]:
Box.test(bond_m2$residuals, lag = 12, type = 'Ljung')
Fitting a time series residuals¶
In [319]:
bond_m3 = arima(aaac, order = c(0,0,1), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [320]:
bond_m3 = arima(aaac, order = c(1,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [323]:
bond_m3 = arima(aaac, order = c(1,0,1), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [ ]:
In [324]:
cat("order =", ar(bond_m2$residuals, method = 'mle')$order, "\n")
In [325]:
bond_m3 = arima(aaac, order = c(12,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [326]:
bond_m3 = arima(aaac, order = c(2,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [327]:
bond_m3 = arima(aaac, order = c(3,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [338]:
bond_m3 = arima(aaac, order = c(4,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [339]:
bond_m3 = arima(aaac, order = c(5,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [340]:
bond_m3 = arima(aaac, order = c(8,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [335]:
bond_m3 = arima(aaac, order = c(2,0,2), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [337]:
bond_m3 = arima(
aaac, order = c(2,0,2),
fixed = c(0, NA, NA, NA, NA),
xreg = bbbc, include.mean = FALSE
)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(bbbc^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Forecast¶
In [332]:
npts = 26
eotr = length(aaac_ts)-npts
h = npts
freq = 52
xreg = as.matrix(bbbc)
aaac_fc_res = plot_auto_arima_forecast_fig(
da_ts=aaac_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
xreg=xreg,
main=NULL, # "Forecasts from ARIMA(2,0,2)\nfor CRSP Equal-Weighted Index",
xlab="Year", ylab="Aaa Bond Yield (using Baa Bond Yield)", # , ylim=c(5, 11)
d = 0,
D = 0,
max.p = 3,
max.q = 4,
max.P = 0,
max.Q = 0,
max.order = 5,
seasonal = FALSE,
method = "ML",
allowmean = FALSE,
stepwise = FALSE,
parallel = TRUE,
num.cores = 2
)
In [333]:
aaac_fc_tb = comb_forecast_res(aaac_fc_res, aaac_ts, eotr, freq)
aaac_fc_tb
In [342]:
par(bg = 'white')
tsdiag(aaac_fc_res$model, gof.lag = 36)
In [ ]:
2-13¶
In [ ]:
In [ ]:
In [ ]:
In [ ]:
2-14¶
In [ ]:
In [ ]:
In [ ]:
In [ ]:
2-15¶
In [ ]:
In [ ]:
In [ ]: